#ifndef SVMTRAIN_CU
#define SVMTRAIN_CU

#include <cuda.h>
#include <cutil.h>
#include <cutil_inline.h>
#include "../common/framework.h"
#include "Cache.h"
//#include "Controller.h"
#include "svmTrain.h"
#include "initialize.h"
#include "firstOrder.h"
#include "secondOrder.h"

void formModel(float* trainingPoints, int nTrainingPoints,
			   int nDimension,
			   float* trainingAlpha,
			   float* trainingLabels,
			   float** supportVectors,
			   int* nSV,
			   float** alpha,
			   float epsilon)
{
	int count = 0;

	for(int i = 0; i < nTrainingPoints; i++) {
		if (trainingAlpha[i] > epsilon) {
			count++;
		}
	}
	*nSV = count;
	printf("%i support vectors found\n", count);
	float* mySupportVectors = *supportVectors;
	mySupportVectors = (float*)malloc(count * nDimension * sizeof(float));
	*supportVectors = mySupportVectors;

	float* myAlpha = *alpha;
	myAlpha = (float*)malloc(count * sizeof(float));
	*alpha = myAlpha;
	int currentOutput = 0;
	for(int i = 0; i < nTrainingPoints; i++) {
		if (trainingAlpha[i] > epsilon) {
			float* sourcePointer = &trainingPoints[i];
			float* destinationPointer = &mySupportVectors[currentOutput];
			for(int j = 0; j < nDimension; j++) {
				*destinationPointer = *sourcePointer;
				sourcePointer += nTrainingPoints;
				destinationPointer += count;
			}
			myAlpha[currentOutput] = trainingAlpha[i] * trainingLabels[i];
			currentOutput++;
		}
	}
}


void performTraining(float* data,
					 int nPoints,
					 int nDimension,
					 float* labels,
					 float** p_alpha,
					 KernelParams* kp,
					 float cost,
					 SelectionHeuristic heuristicMethod,
					 float epsilon,
					 float tolerance,
					 float* transposedData,
					 int ellWidth,
					 float* ellEntry,
					 int* ellColumn,
					 int* ellRowLen)
{
	float cEpsilon = cost - epsilon;
	//TODO 2. remove progress CONTROLLER here
	//Controller progress(2.0, heuristicMethod, 64, nPoints);

	float parameterA;
	float parameterB;
	float parameterC;
	if (kp->kType == GAUSSIAN) {
		parameterA = -kp->gamma;
	} else if (kp->kType == POLYNOMIAL) {
		parameterA = kp->gamma;
		parameterB = kp->coef0;
		parameterC = kp->degree;
	} else if (kp->kType == SIGMOID) {
		parameterA = kp->gamma;
		parameterB = kp->coef0;
	}
	printf("> Cost: %f, Tolerance: %f, Epsilon: %f\n", cost, tolerance, epsilon);

	float* devData;
	float* devTransposedData;
	size_t devDataPitch;
	size_t devTransposedDataPitch;

	CUDA_SAFE_CALL(cudaMallocPitch((void**)&devData, &devDataPitch, nPoints*sizeof(float)/*in bytes*/, nDimension));
	CUDA_SAFE_CALL(cudaMemcpy2D(devData, devDataPitch, data, nPoints*sizeof(float)/*pitch of source mem*/, nPoints*sizeof(float)/*width in bytes*/, nDimension ,cudaMemcpyHostToDevice));
	
	CUDA_SAFE_CALL(cudaMallocPitch((void**)&devTransposedData, &devTransposedDataPitch, nDimension*sizeof(float), nPoints));
	CUDA_SAFE_CALL(cudaMemcpy2D(devTransposedData, devTransposedDataPitch, transposedData, nDimension*sizeof(float), nDimension*sizeof(float), nPoints, cudaMemcpyHostToDevice));

	float* alpha = (float*)malloc(sizeof(float) * nPoints);
	*p_alpha = alpha;

	float* devLabels;
	CUDA_SAFE_CALL(cudaMalloc((void**)&devLabels, nPoints*sizeof(float)));
	CUDA_SAFE_CALL(cudaMemcpy(devLabels, labels, nPoints*sizeof(float), cudaMemcpyHostToDevice));

	float* devKernelDiag;
	CUDA_SAFE_CALL(cudaMalloc((void**)&devKernelDiag, nPoints*sizeof(float)));

	float* devAlpha;
	CUDA_SAFE_CALL(cudaMalloc((void**)&devAlpha, nPoints*sizeof(float)));

	float* devF;
	CUDA_SAFE_CALL(cudaMalloc((void**)&devF, nPoints*sizeof(float)));

	void* devResult;
	CUDA_SAFE_CALL(cudaMalloc(&devResult, 8*sizeof(float)));
	float* hostResult = (float*)malloc(8*sizeof(float));

	int blockWidth = (nPoints+(BLOCKSIZE-1))/BLOCKSIZE;

	float* devLocalFsRL;
	CUDA_SAFE_CALL(cudaMalloc((void**)&devLocalFsRL, blockWidth*sizeof(float)));
	float* devLocalFsRH;
	CUDA_SAFE_CALL(cudaMalloc((void**)&devLocalFsRH, blockWidth*sizeof(float)));
	int* devLocalIndicesRL;
	CUDA_SAFE_CALL(cudaMalloc((void**)&devLocalIndicesRL, blockWidth*sizeof(int)));
	int* devLocalIndicesRH;
	CUDA_SAFE_CALL(cudaMalloc((void**)&devLocalIndicesRH, blockWidth*sizeof(int)));

	float* devLocalObjsMaxObj;
	CUDA_SAFE_CALL(cudaMalloc((void**)&devLocalObjsMaxObj, blockWidth*sizeof(float)));
	int* devLocalIndicesMaxObj;
	CUDA_SAFE_CALL(cudaMalloc((void**)&devLocalIndicesMaxObj, blockWidth*sizeof(int)));

	//ELLPACK data allocation
	float* devEllEntry;
	int* devEllColumn;
	int* devEllRowLen;
	size_t devEllEntryPitch,devEllColPitch;
	CUDA_SAFE_CALL(cudaMallocPitch((void**) &devEllEntry, &devEllEntryPitch, ellWidth*sizeof(float), nPoints));
	CUDA_SAFE_CALL(cudaMemcpy2D(devEllEntry, devEllEntryPitch, ellEntry, ellWidth*sizeof(float), ellWidth*sizeof(float), nPoints, cudaMemcpyHostToDevice));
	cudaChannelFormatDesc channelDescEllEntry = cudaCreateChannelDesc<float>();
	ellEntryTex.normalized = 0;
	ellEntryTex.filterMode = cudaFilterModePoint;
	ellEntryTex.addressMode[0] = cudaAddressModeClamp;
	ellEntryTex.addressMode[1] = cudaAddressModeClamp;
	CUDA_SAFE_CALL(cudaBindTexture2D(NULL, ellEntryTex, devEllEntry, channelDescEllEntry, ellWidth, nPoints, devEllEntryPitch));

	CUDA_SAFE_CALL(cudaMallocPitch((void**) &devEllColumn, &devEllColPitch, ellWidth*sizeof(int), nPoints));
	CUDA_SAFE_CALL(cudaMemcpy2D(devEllColumn, devEllColPitch, ellColumn, ellWidth*sizeof(int), ellWidth*sizeof(int), nPoints, cudaMemcpyHostToDevice));
	cudaChannelFormatDesc channelDescEllColumn = cudaCreateChannelDesc<int>();
	ellColumnTex.normalized = 0;
	ellColumnTex.filterMode = cudaFilterModePoint;
	ellColumnTex.addressMode[0] = cudaAddressModeClamp;
	ellColumnTex.addressMode[1] = cudaAddressModeClamp;
	CUDA_SAFE_CALL(cudaBindTexture2D(NULL, ellColumnTex, devEllColumn, channelDescEllColumn, ellWidth, nPoints, devEllColPitch));

	CUDA_SAFE_CALL(cudaMalloc((void**) &devEllRowLen, nPoints*sizeof(int)));
	CUDA_SAFE_CALL(cudaMemcpy(devEllRowLen, ellRowLen, nPoints*sizeof(int), cudaMemcpyHostToDevice));
	//cudaChannelFormatDesc channelDescEllRowLen = cudaCreateChannelDesc<int>();
	//ellRowLenTex.normalized = 0;
	//ellRowLenTex.filterMode = cudaFilterModePoint;
	//ellRowLenTex.addressMode[0] = cudaAddressModeClamp;
	//CUDA_SAFE_CALL(cudaBindTexture(NULL, ellRowLenTex, devEllRowLen, channelDescEllRowLen, nPoints*sizeof(int)));
	CUDA_SAFE_CALL(cudaBindTexture(NULL, ellRowLenTex, devEllRowLen, nPoints*sizeof(int)));

	unsigned int remainingMemory;
	unsigned int totalMemory;
	cuMemGetInfo(&remainingMemory, &totalMemory);

	int sizeOfCacheInRow = remainingMemory/((int)devDataPitch);
	//TODO 4. If I try to grab all the memory available, it'll fail
	sizeOfCacheInRow = (int)((float)sizeOfCacheInRow*0.95);
	if (nPoints < sizeOfCacheInRow) {
		sizeOfCacheInRow = nPoints;
	}

#ifdef __DEVICE_EMULATION__
	sizeOfCache = nPoints;
#endif

	printf("%d bytes of memory found on device, %d bytes currently free\n", totalMemory, remainingMemory);
	printf("%d rows of kernel matrix will be cached (%d bytes per row)\n", sizeOfCacheInRow, (int)devDataPitch);

	float* devCache;
	size_t cachePitch;
	CUDA_SAFE_CALL(cudaMallocPitch((void**)&devCache, &cachePitch, nPoints*sizeof(float), sizeOfCacheInRow));
	Cache kernelCache(nPoints, sizeOfCacheInRow);
	int devCachePitchInFloats = (int)cachePitch/(sizeof(float));

	CUERR
	printf("Allocated arrays on GPU\n");

	dim3 threadsLinear(BLOCKSIZE);
	dim3 blocksLinear(blockWidth);

	int devDataPitchInFloats = ((int)devDataPitch)/(sizeof(float));
	int devTransposedDataPitchInFloats = ((int)devTransposedDataPitch)/(sizeof(float));
	
	//TODO 6. initialize alphaI, Fi, and diagonal Kernal pre-comute value(Kii)
	launchInitialization(devData, devDataPitchInFloats, nPoints, nDimension, kp->kType, parameterA, parameterB, parameterC, devKernelDiag, devAlpha, devF, devLabels, blocksLinear, threadsLinear);
	
	CUERR
	printf("Initialization complete\n");

	//Choose initial points
	float bLow = 1;
	float bHigh = -1;
	int iteration = 0;
	int iLow = -1;
	int iHigh = -1;
	for (int i = 0; i < nPoints; i++) {
		if (labels[i] < 0) {
			if (iLow == -1) {
				iLow = i;
				if (iHigh != -1) break;
			}
		} else {
			if (iHigh == -1) {
				iHigh = i;
				if (iLow != -1) break;
			}
		}
	}

	dim3 singletonThreads(1);
	dim3 singletonBlocks(1);
	launchTakeFirstStep(devResult, devKernelDiag, devData, devDataPitchInFloats, devAlpha, cost, nDimension, iLow, iHigh, kp->kType, parameterA, parameterB, parameterC, singletonBlocks, singletonThreads);
	CUDA_SAFE_CALL(cudaMemcpy((void*)hostResult, devResult, 8*sizeof(float), cudaMemcpyDeviceToHost));

	float alpha2Old = *(hostResult + 0);
	float alpha1Old = *(hostResult + 1);
	bLow = *(hostResult + 2);
	bHigh = *(hostResult + 3);
	float alpha2New = *(hostResult + 6);
	float alpha1New = *(hostResult + 7);

	float alpha1Diff = alpha1New - alpha1Old;
	float alpha2Diff = alpha2New - alpha2Old;

	int iLowCacheIndex;
	int iHighCacheIndex;
	bool iLowCompute;
	bool iHighCompute; 

	dim3 reduceThreads(BLOCKSIZE);

	for (iteration = 1 ;; iteration++) {
		if (bLow <= bHigh + 2*tolerance) {
			break; //Convergence!!
		}

		if ((iteration & 0x7ff) == 0) {
			printf("iteration: %d; gap: %f; 2*tol: %f\n",iteration, bLow - bHigh, 2.0f*tolerance);
		}
		//TODO 2. remove progress CONTROLLER here
		//if ((iteration & 0x7f) == 0) {
		//	heuristicMethod = progress.getMethod();
		//}

		kernelCache.findData(iHigh, iHighCacheIndex, iHighCompute);
		kernelCache.findData(iLow, iLowCacheIndex, iLowCompute);

		if (heuristicMethod == FIRSTORDER) {
			launchFirstOrder(iLowCompute, iHighCompute, kp->kType, nPoints, nDimension, blocksLinear, threadsLinear, reduceThreads, devData, devDataPitchInFloats, devTransposedData, devTransposedDataPitchInFloats, devLabels, epsilon, cEpsilon, devAlpha, devF, alpha1Diff * labels[iHigh], alpha2Diff * labels[iLow], iLow, iHigh, parameterA, parameterB, parameterC, devCache, devCachePitchInFloats, iLowCacheIndex, iHighCacheIndex, devLocalIndicesRL, devLocalIndicesRH, devLocalFsRH, devLocalFsRL, devKernelDiag, devResult, cost);
		} else if(heuristicMethod == SECONDORDER){
			launchSecondOrder(iLowCompute, iHighCompute, kp->kType, nPoints, nDimension, blocksLinear, threadsLinear, reduceThreads, devData, devDataPitchInFloats, devTransposedData, devTransposedDataPitchInFloats, devLabels, epsilon, cEpsilon, devAlpha, devF, alpha1Diff * labels[iHigh], alpha2Diff * labels[iLow], iLow, iHigh, parameterA, parameterB, parameterC, &kernelCache, devCache, devCachePitchInFloats, iLowCacheIndex, iHighCacheIndex, devLocalIndicesRH, devLocalFsRH, devLocalFsRL, devLocalIndicesMaxObj, devLocalObjsMaxObj, devKernelDiag, devResult, hostResult, cost, iteration);
		} else {
			printf(">!!!Wrong heuristic method.\n");
			exit(1);
		}
		//printf("iteration: %i\n", iteration);
		CUDA_SAFE_CALL(cudaMemcpy((void*)hostResult, devResult, 8*sizeof(float), cudaMemcpyDeviceToHost));

		alpha2Old = *(hostResult + 0);
		alpha1Old = *(hostResult + 1);
		bLow = *(hostResult + 2);
		bHigh = *(hostResult + 3);
		alpha2New = *(hostResult + 4);
		alpha1New = *(hostResult + 5);
		iLow = *((int*)hostResult + 6);
		iHigh = *((int*)hostResult + 7);
		
		alpha1Diff = alpha1New - alpha1Old;
		alpha2Diff = alpha2New - alpha2Old;
		/*TODO 2. remove progress CONTROLLER here*/
		//progress.addIteration(bLow-bHigh);
	}

	printf("%d iterations\n", iteration);
	printf("bLow: %f, bHigh: %f\n", bLow, bHigh);
	kp->b = (bLow + bHigh) / 2;
	kernelCache.printStatistics();
	//kernelCache.printCache();
	CUDA_SAFE_CALL(cudaMemcpy((void*)alpha, devAlpha, nPoints*sizeof(float), cudaMemcpyDeviceToHost));
	cudaUnbindTexture(ellEntryTex);
	cudaFree(devEllEntry);
	cudaUnbindTexture(ellColumnTex);
	cudaFree(devEllColumn);
	cudaUnbindTexture(ellRowLenTex);
	cudaFree(devEllRowLen);
	cudaFree(devData);
	cudaFree(devTransposedData);
	cudaFree(devLabels);
	cudaFree(devAlpha);
	cudaFree(devF);
	cudaFree(devCache);
	cudaFree(devLocalIndicesRL);
	cudaFree(devLocalIndicesRH);
	cudaFree(devLocalFsRH);
	cudaFree(devLocalFsRL);
	cudaFree(devKernelDiag);
	cudaFree(devResult);
	cudaFree(devLocalIndicesMaxObj);
	cudaFree(devLocalObjsMaxObj);
	free(hostResult);
	system("PAUSE");
}
#endif